using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// Calculates an approximation result to a given definite integral I = integral of F over (A,B).
    /// </summary>
    [Serializable]
    public static class DqngClass
    {
        /// <summary>
        /// Calculates an approximation result to a given definite integral I = integral of F over (A,B).
        /// </summary>
        /// <param name="f">Function subprogram defining the integrand function.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="epsabs">Absolute accuracy requested.</param>
        /// <param name="epsrel">Relative accuracy requested.</param>
        /// <param name="abserr">Absolute accuracy requested.</param>
        /// <param name="neval">The needed function evaluations.</param>
        /// <param name="ier">The error code.</param>
        /// <returns>Approximation to the integral.</returns>
        public static double Dqng(HardRealFunction f, double a, double b, double epsabs, double epsrel,
                                  ref double abserr, ref int neval, ref int ier)
        {
            double[] x1 = {
                              0.97390652851717172008, 0.86506336668898451073, 0.67940956829902440623,
                              0.43339539412924719080, 0.14887433898163121088
                          };
            double[] w10 = {
                               0.06667134430868813759, 0.14945134915058059315, 0.21908636251598204400,
                               0.26926671930999635509, 0.29552422471475287017
                           };
            double[] x2 = {
                              0.99565716302580808074, 0.93015749135570822600, 0.78081772658641689706,
                              0.56275713466860468334, 0.29439286270146019813
                          };
            double[] w21a = {
                                0.03255816230796472748, 0.07503967481091995277, 0.10938715880229764190,
                                0.13470921731147332593, 0.14773910490133849137
                            };
            double[] w21b = {
                                0.01169463886737187428, 0.05475589657435199603, 0.09312545458369760554,
                                0.12349197626206585108, 0.14277593857706008080, 0.14944555400291690566
                            };
            double[] x3 = {
                              0.99933336090193208139, 0.98743340290808886980, 0.95480793481426629926,
                              0.90014869574832829363, 0.82519831498311415085, 0.73214838898930498261,
                              0.62284797053772523864, 0.49947957407105649995, 0.36490166134658076804,
                              0.22225491977660129650, 0.07465061746138332204
                          };
            double[] w43a = {
                                0.01629673428966656492, 0.03752287612086950146, 0.05469490205825544215,
                                0.06735541460947808608, 0.07387019963239395343, 0.00576855605976979618,
                                0.02737189059324884208, 0.04656082691042883074, 0.06174499520144256450,
                                0.07138726726869339777
                            };
            double[] w43b = {
                                0.00184447764021241410, 0.01079868958589165174, 0.02189536386779542810,
                                0.03259746397534568944, 0.04216313793519181185, 0.05074193960018457778,
                                0.05837939554261924838, 0.06474640495144588554, 0.06956619791235648453,
                                0.07282444147183320815, 0.07450775101417511827, 0.07472214751740300559
                            };
            double[] x4 = {
                              0.99990297726272923449, 0.99798989598667874543, 0.99217549786068722281,
                              0.98135816357271277357, 0.96505762385838461913, 0.94316761313367059682,
                              0.91580641468550720959, 0.88322165777131650137, 0.84571074846241566661,
                              0.80355765803523098279, 0.75700573068549555833, 0.70627320978732181982,
                              0.65158946650117792253, 0.59322337405796108888, 0.53149360597083193229,
                              0.46676362304202284487, 0.39942484785921880473, 0.32987487710618828827,
                              0.25850355920216155180, 0.18569539656834665202, 0.11184221317990746817,
                              0.03735212339461987081
                          };
            double[] w87a = {
                                0.00814837738414917290, 0.01876143820156282224, 0.02734745105005228616,
                                0.03367770731163793005, 0.03693509982042790761, 0.00288487243021153050,
                                0.01368594602271270189, 0.02328041350288831112, 0.03087249761171335868,
                                0.03569363363941877072, 0.00091528334520224136, 0.00539928021930047137,
                                0.01094767960111893113, 0.01629873169678733526, 0.02108156888920383511,
                                0.02537096976925382724, 0.02918969775647575250, 0.03237320246720278969,
                                0.03478309895036514275, 0.03641222073135178756, 0.03725387550304770854
                            };
            double[] w87b = {
                                0.00027414556376207235, 0.00180712415505794295, 0.00409686928275916486,
                                0.00675829005184737870, 0.00954995767220164654, 0.01232944765224485369,
                                0.01501044734638895238, 0.01754896798624319110, 0.01993803778644088820,
                                0.02219493596101228680, 0.02433914712600080547, 0.02637450541483920724,
                                0.02828691078877120066, 0.03005258112809269532, 0.03164675137143992940,
                                0.03305041341997850329, 0.03425509970422606179, 0.03526241266015668103,
                                0.03607698962288870119, 0.03669860449845609450, 0.03712054926983257611,
                                0.03733422875193504032, 0.03736107376267902341
                            };

            double[] fv1 = new double[5];
            double[] fv2 = new double[5];
            double[] fv3 = new double[5];
            double[] fv4 = new double[5];
            double[] savfun = new double[21];
            double absc;
            double centr;
            double dhlgth;
            double fcentr;
            double fval;
            double fval1;
            double fval2;
            double hlgth;
            double result;
            double res10;
            double res21 = 0;
            double res43 = 0;
            double res87;
            double resabs = 0;
            double resasc = 0;
            double reskh;
            int ipx = 0;
            int k;
            int l;

            result = 0.0;
            abserr = 0.0;
            neval = 0;
            ier = 6;
            if ((epsabs < 0.0) && (epsrel < 0.0))
            {
                return result;
            }
            hlgth = 0.5 * (b - a);
            dhlgth = Math.Abs(hlgth);
            centr = 0.5 * (a + b);
            fcentr = f.SolveAt(centr);
            neval = 21;
            ier = 1;

            for (l = 1; l <= 3; l++)
            {
                switch (l)
                {
                    case 1:
                        res10 = 0.0;
                        res21 = w21b[5] * fcentr;
                        resabs = w21b[5] * Math.Abs(fcentr);
                        for (k = 0; k < 5; k++)
                        {
                            absc = hlgth * x1[k];
                            fval1 = f.SolveAt(centr + absc);
                            fval2 = f.SolveAt(centr - absc);
                            fval = fval1 + fval2;
                            res10 += (w10[k] * fval);
                            res21 += (w21a[k] * fval);
                            resabs += (w21a[k] * (Math.Abs(fval1) + Math.Abs(fval2)));
                            savfun[k] = fval;
                            fv1[k] = fval1;
                            fv2[k] = fval2;
                        }
                        ipx = 4;
                        for (k = 0; k < 5; k++)
                        {
                            ipx++;
                            absc = hlgth * x2[k];
                            fval1 = f.SolveAt(centr + absc);
                            fval2 = f.SolveAt(centr - absc);
                            fval = fval1 + fval2;
                            res21 += (w21b[k] * fval);
                            resabs += (w21b[k] * (Math.Abs(fval1) + Math.Abs(fval2)));
                            savfun[ipx] = fval;
                            fv3[k] = fval1;
                            fv4[k] = fval2;
                        }
                        result = res21 * hlgth;
                        resabs *= dhlgth;
                        reskh = 0.5 * res21;
                        resasc = w21b[5] * Math.Abs(fcentr - reskh);
                        for (k = 0; k < 5; k++)
                        {
                            resasc += (w21a[k] * (Math.Abs(fv1[k] - reskh) + Math.Abs(fv2[k] - reskh)) +
                                       w21b[k] * (Math.Abs(fv3[k] - reskh) + Math.Abs(fv4[k] - reskh)));
                        }
                        abserr = Math.Abs((res21 - res10) * hlgth);
                        resasc *= dhlgth;
                        break;
                    case 2:
                        res43 = w43b[11] * fcentr;
                        neval = 43;
                        for (k = 0; k < 10; k++)
                        {
                            res43 += (savfun[k] * w43a[k]);
                        }
                        for (k = 0; k < 11; k++)
                        {
                            ipx++;
                            absc = hlgth * x3[k];
                            fval = f.SolveAt(centr + absc) + f.SolveAt(centr - absc);
                            res43 += (fval * w43b[k]);
                            savfun[ipx] = fval;
                        }
                        result = res43 * hlgth;
                        abserr = Math.Abs((res43 - res21) * hlgth);
                        break;
                    case 3:
                        res87 = w87b[22] * fcentr;
                        neval = 87;
                        for (k = 0; k < 21; k++)
                        {
                            res87 += (savfun[k] * w87a[k]);
                        }
                        for (k = 0; k < 22; k++)
                        {
                            absc = hlgth * x4[k];
                            res87 += w87b[k] * (f.SolveAt(centr + absc) + f.SolveAt(centr - absc));
                        }
                        result = res87 * hlgth;
                        abserr = Math.Abs((res87 - res43) * hlgth);
                        break;
                }
                if ((resasc != 0.0) && (abserr != 0.0))
                {
                    abserr = resasc *
                             (((1.0) < (Math.Pow(200.0 * (abserr) / resasc, 1.5)))
                                  ? (1.0)
                                  : (Math.Pow(200.0 * (abserr) / resasc, 1.5)));
                }
                if (resabs > Double.MinValue / (50.0 * Double.Epsilon))
                {
                    abserr = ((((Double.Epsilon * 50.0) * resabs) > (abserr))
                                  ? ((Double.Epsilon * 50.0) * resabs)
                                  : (abserr));
                }
                if (abserr <= (((epsabs) > (epsrel * Math.Abs(result))) ? (epsabs) : (epsrel * Math.Abs(result))))
                {
                    ier = 0;
                }
                if (ier == 0)
                {
                    break;
                }
            }
            return result;
        }
    }
}